BARRETO et al. 2022: Habitat loss and biodiversity gain (BIOCON-S-22-00746)

Model formulas

For each of our response variables we run three candidate models for each scale:

Model 0 <- ~1, absence of effect

Model 1 <- percentage of forest cover at 3 and 5km radius

Model1.2 <- non-linear forest cover at 3 and 5km radius (quadratic model)

Loading data and gathering into one datasets calculated in script ‘data_wrang_GIT.Rmd’, content: 1) “beetles”: raw data for species occurrence per landscape/site; 2) “hab.data”: data containing classes of habitat association; 3) “env.data”: environmental variables, forest cover percent at 3 and 5km radius buffer 4) “div_hab”: diversity (alpha, beta and gamma) between classification habitat association 5) “rc.df”: Beta RC calculations 6) “ab_hab”: mean and median abundance between groups of habitat association 7) “st.hab_ab”: abundance data between classification habitat association per site

Landscape gamma_FS gamma_NFS alpha_FS alpha_NFS betad_FS betad_NFS ab.land_FS ab.land_NFS brc.fsNA brc.nfsNA brc.abund.fsNA brc.abund.nfsNA fc_3km fc_5km
148 11 9 3.75 2.25 7.25 6.75 66 40 -0.2989286 0.1920000 -0.0067210 0.4313647 38.52523 43.86337
215 19 7 5.62 2.00 13.38 5.00 174 51 0.0989286 -0.3042857 0.2793534 0.4023071 30.15485 25.23160
263 16 13 6.00 4.50 10.00 8.50 393 214 -0.5085714 -0.0057143 0.6216216 0.4039754 10.10641 13.41180
266 11 8 3.25 1.38 7.75 6.62 70 12 -0.2342857 0.1352381 0.1091091 0.1504362 30.87525 31.16090
282 18 10 6.62 4.00 11.38 6.00 220 211 -0.2835714 -0.1852381 0.2157157 0.6922637 15.26270 16.35308
291 9 3 3.71 1.71 5.29 1.29 150 49 -0.3938095 -0.5809524 0.0873731 -0.4176081 45.74294 45.13213
317 12 4 5.00 1.38 7.00 2.62 152 56 -0.4033333 -0.5364286 0.1798584 -0.4283212 48.75805 47.85138
329 21 13 8.75 5.50 12.25 7.50 373 258 -0.0692593 0.1207407 0.9381882 0.4310382 14.73791 13.62588
333 15 14 4.62 4.50 10.38 9.50 159 218 -0.1096429 0.1907143 0.6733877 0.6187259 13.71403 14.50960
335 15 8 5.75 2.38 9.25 5.62 140 113 -0.2660714 -0.3266667 0.3353353 -0.9401783 23.87913 19.21670
359 17 16 6.25 5.25 10.75 10.75 152 130 -0.1721429 0.3191667 0.7743815 0.8145288 18.10208 17.62431
399 21 13 8.00 5.00 13.00 8.00 433 168 -0.0852381 0.2525926 0.9580056 0.7037752 29.73616 24.97851

RESULTS

We found 28 species to be forest specialists, those that come from the Atlantic Forest domain and exclusively collected in forested habitats; while 23 were non-forest specialists, species from a broader domain and/or less restrictive to forest habitats (hereafter FS and NFS, respectively).

FOREST SPECIALISTS

Gamma diversity

Modelling

g0 <- glm(gamma_FS ~ 1, data=data_div, family= Gamma (link = "identity"))
g5k <- glm(gamma_FS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))
g5kq <- glm(gamma_FS ~ fc_5km +  I(fc_5km^2), data=data_div, family= Gamma (link = "identity"))
g3k <- glm(gamma_FS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))
g3kq <- glm(gamma_FS ~ fc_3km +  I(fc_3km^2), data=data_div, family= Gamma (link = "identity"))

## AIC-based model selection
(g.fs.aic <- AICctab(g0, g5k, g5kq, g3k, g3kq,
                      base=T,weights=T, logLik=T))
##      logLik AICc  dLogLik dAICc df weight
## g5k  -28.2   65.4   5.1     0.0 3  0.751 
## g3k  -30.0   68.9   3.3     3.5 3  0.127 
## g5kq -28.1   69.9   5.2     4.6 4  0.076 
## g0   -33.3   71.9   0.0     6.6 2  0.028 
## g3kq -29.6   72.9   3.7     7.6 4  0.017

Summary top model

summary(g5k)
## 
## Call:
## glm(formula = gamma_FS ~ fc_5km, family = Gamma(link = "identity"), 
##     data = data_div)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.24788  -0.14644  -0.02457   0.14244   0.30719  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 21.46708    2.02439  10.604 9.26e-07 ***
## fc_5km      -0.23174    0.05885  -3.938  0.00278 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.03547038)
## 
##     Null deviance: 0.80342  on 11  degrees of freedom
## Residual deviance: 0.34475  on 10  degrees of freedom
## AIC: 62.361
## 
## Number of Fisher Scoring iterations: 3
g.fs <- g5k# rename top model to save and plot below

Residuals

sim <- simulateResiduals(g5k)
plot(sim)

testDispersion(sim)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 1.1271, p-value = 0.672
## alternative hypothesis: two.sided
hnp::hnp(g5k)
## Gamma model

boot::glm.diag.plots(g5k) ## for glm poisson or neg binom

Spatial autocorrelation All Moran’s I tests with our response variables to verify spatial auto-correlation resulted in non-significant which means non spatially auto-correlated.

testSpatialAutocorrelation(sim, x =  env.data$x, y = env.data$y)

## 
##  DHARMa Moran's I test for distance-based autocorrelation
## 
## data:  sim
## observed = 0.013940, expected = -0.090909, sd = 0.057965, p-value =
## 0.07048
## alternative hypothesis: Distance-based autocorrelation

Prediction

Alpha diversity

Modelling

a0 <- glm(alpha_FS ~ 1, data=data_div, family= Gamma (link = "identity"))
a5k <- glm(alpha_FS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))
a5kq <- glm(alpha_FS ~ fc_5km +  I(fc_5km^2), data=data_div, family= Gamma (link = "identity"))
a3k <- glm(alpha_FS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))
a3kq <- glm(alpha_FS ~ fc_3km +  I(fc_3km^2), data=data_div, family= Gamma (link = "identity"))


## AIC-based model selection
(a.fs.aic <- AICctab(a0, a5k, a5kq,
                     a3k, a3kq, 
                     base=T,weights=T, logLik=T))
##      logLik AICc  dLogLik dAICc df weight
## a5k  -19.5   48.0   3.0     0.0 3  0.538 
## a3k  -20.5   50.0   2.0     2.0 3  0.200 
## a0   -22.4   50.2   0.0     2.2 2  0.175 
## a5kq -19.2   52.1   3.2     4.1 4  0.068 
## a3kq -20.5   54.7   2.0     6.7 4  0.019

Summary top model

summary(a5k)
## 
## Call:
## glm(formula = alpha_FS ~ fc_5km, family = Gamma(link = "identity"), 
##     data = data_div)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.43968  -0.11810  -0.03523   0.09321   0.36169  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  7.56549    0.96276   7.858 1.38e-05 ***
## fc_5km      -0.07518    0.02853  -2.635    0.025 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.06116683)
## 
##     Null deviance: 1.00974  on 11  degrees of freedom
## Residual deviance: 0.62036  on 10  degrees of freedom
## AIC: 44.989
## 
## Number of Fisher Scoring iterations: 4
a.fs <- a5k# rename top model to save and plot below

Residuals

sim <- simulateResiduals(a5k)
plot(sim)

testDispersion(sim)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 1.1023, p-value = 0.672
## alternative hypothesis: two.sided
hnp::hnp(a5k)
## Gamma model

boot::glm.diag.plots(a5k)

Spatial autocorrelation

testSpatialAutocorrelation(sim, x =  env.data$x, y = env.data$y)

## 
##  DHARMa Moran's I test for distance-based autocorrelation
## 
## data:  sim
## observed = 0.0054053, expected = -0.0909091, sd = 0.0566974, p-value =
## 0.08937
## alternative hypothesis: Distance-based autocorrelation

Prediction

Beta diversity additive

Modelling

### LM Beta for regular model selection
b0 <- glm(betad_FS~ 1, data=data_div, family= Gamma (link = "identity"))
b5k <- glm(betad_FS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))
b5kq <- glm(betad_FS ~ fc_5km +  I(fc_5km^2), data=data_div, family= Gamma (link = "identity"))
b3k <- glm(betad_FS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))
b3kq <- glm(betad_FS ~ fc_3km +  I(fc_3km^2), data=data_div, family= Gamma (link = "identity"))

## AIC-based model selection
(b.fs.aic <- AICctab(b0, b5k, b5kq, 
                     b3k, b3kq, 
                     base=T,weights=T, logLik=T))
##      logLik AICc  dLogLik dAICc df weight
## b5k  -22.5   54.0   5.6     0.0 3  0.710 
## b5kq -21.8   57.3   6.3     3.3 4  0.135 
## b3k  -24.4   57.7   3.7     3.8 3  0.108 
## b3kq -23.3   60.2   4.8     6.2 4  0.031 
## b0   -28.1   61.5   0.0     7.6 2  0.016

Summary top model

summary(b5k)
## 
## Call:
## glm(formula = betad_FS ~ fc_5km, family = Gamma(link = "identity"), 
##     data = data_div)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.24292  -0.15104  -0.01746   0.05238   0.31092  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 13.93904    1.28413  10.855 7.46e-07 ***
## fc_5km      -0.15787    0.03686  -4.283   0.0016 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.03504058)
## 
##     Null deviance: 0.83647  on 11  degrees of freedom
## Residual deviance: 0.33078  on 10  degrees of freedom
## AIC: 50.972
## 
## Number of Fisher Scoring iterations: 4
b.fs <- b5k # rename top model to save and plot below

Residuals

sim <- simulateResiduals(b5k)
plot(sim)

testDispersion(sim) # ok

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 1.1273, p-value = 0.624
## alternative hypothesis: two.sided
hnp::hnp(b5k) ## for glm Gamma or neg binom (maybe also others)
## Gamma model

boot::glm.diag.plots(b5k) ## for glm Gamma or neg binom

Spatial autocorrelation

testSpatialAutocorrelation(sim, x =  env.data$x, y = env.data$y)

## 
##  DHARMa Moran's I test for distance-based autocorrelation
## 
## data:  sim
## observed = 0.021701, expected = -0.090909, sd = 0.057512, p-value =
## 0.05023
## alternative hypothesis: Distance-based autocorrelation

Prediction

bRaup-Crick (presence/absence-based)

Modelling

### LM gamma for regular model selection
rc0 <- lm(brc.fsNA ~ 1, data=data_div)
rc5k <- lm(brc.fsNA ~ fc_5km, data=data_div)
rc5kq <- lm(brc.fsNA ~ fc_5km +  I(fc_5km^2), data=data_div)
rc3k <- lm(brc.fsNA ~ fc_3km, data=data_div)
rc3kq <- lm(brc.fsNA ~ fc_3km +  I(fc_3km^2), data=data_div)

## AIC-based model selection
(rc.fs.aic <- AICctab(rc0, rc5k, rc5kq, 
                      rc3k, rc3kq,
                      base=T,weights=T, logLik=T))
##       logLik AICc dLogLik dAICc df weight
## rc0    4.8   -4.2  0.0     0.0  2  0.419 
## rc3kq  8.4   -3.0  3.6     1.2  4  0.227 
## rc5k   5.5   -2.0  0.7     2.2  3  0.139 
## rc5kq  7.8   -1.9  3.0     2.4  4  0.129 
## rc3k   5.0   -1.1  0.2     3.2  3  0.086

Summary of 2nd best model, as the top is of absence of effect (~1)

summary(rc3kq)
## 
## Call:
## lm(formula = brc.fsNA ~ fc_3km + I(fc_3km^2), data = data_div)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.16680 -0.10052 -0.01060  0.07993  0.21297 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -0.6934193  0.2287612  -3.031   0.0142 *
## fc_3km       0.0426487  0.0178385   2.391   0.0405 *
## I(fc_3km^2) -0.0007772  0.0003009  -2.583   0.0296 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1391 on 9 degrees of freedom
## Multiple R-squared:  0.449,  Adjusted R-squared:  0.3266 
## F-statistic: 3.667 on 2 and 9 DF,  p-value: 0.06842
rc1.fs <- rc0# rename top model to save and plot below

Residuals

sim <- simulateResiduals(rc3kq)
plot(sim)

testDispersion(sim) # ok

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.82544, p-value = 0.816
## alternative hypothesis: two.sided
hnp::hnp(rc3kq)
## Gaussian model (lm object)

Spatial autocorrelation

testSpatialAutocorrelation(sim, x =  env.data$x, y = env.data$y)

## 
##  DHARMa Moran's I test for distance-based autocorrelation
## 
## data:  sim
## observed = -0.149708, expected = -0.090909, sd = 0.057697, p-value =
## 0.3082
## alternative hypothesis: Distance-based autocorrelation

Prediction

Total abundance

Modelling

### LM "gamma"  abundanec for regular model selection
gab0 <- glm.nb(ab.land_FS ~ 1, data=data_div)
gab5k <- glm.nb(ab.land_FS ~ fc_5km, data=data_div)
gab5kq <- glm.nb(ab.land_FS ~ fc_5km +  I(fc_5km^2), data=data_div)
gab3k <- glm.nb(ab.land_FS ~ fc_3km, data=data_div)
gab3kq <- glm.nb(ab.land_FS ~ fc_3km +  I(fc_3km^2), data=data_div)

## AIC-based model selection
(gab.fs.aic <- AICctab(gab0, gab5k, gab5kq, 
                       gab3k, gab3kq,
                       base=T,weights=T, logLik=T))
##        logLik AICc  dLogLik dAICc df weight
## gab5k  -70.6  150.1   2.1     0.0 3  0.404 
## gab0   -72.7  150.7   0.0     0.5 2  0.311 
## gab3k  -71.2  151.4   1.5     1.2 3  0.218 
## gab5kq -70.4  154.6   2.2     4.5 4  0.044 
## gab3kq -71.1  155.8   1.6     5.7 4  0.023

Summary 2nd best model, as top model was the one of absence of effect (~1):

summary(gab5k)
## 
## Call:
## glm.nb(formula = ab.land_FS ~ fc_5km, data = data_div, init.theta = 4.550032026, 
##     link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6556  -0.9899  -0.2919   0.6493   1.8325  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  5.93867    0.31886   18.62   <2e-16 ***
## fc_5km      -0.02505    0.01108   -2.26   0.0238 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(4.55) family taken to be 1)
## 
##     Null deviance: 17.357  on 11  degrees of freedom
## Residual deviance: 12.410  on 10  degrees of freedom
## AIC: 147.13
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  4.55 
##           Std. Err.:  1.84 
## 
##  2 x log-likelihood:  -141.134
gab.fs <- gab5k# rename top model to save and plot below

Residuals

sim <- simulateResiduals(gab5k)
plot(sim)

testDispersion(sim)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 1.1526, p-value = 0.496
## alternative hypothesis: two.sided
hnp::hnp(gab5k)
## Negative binomial model (using MASS package)

boot::glm.diag.plots(gab5k)

Spatial autocorrelation

testSpatialAutocorrelation(sim, x =  env.data$x, y = env.data$y)

## 
##  DHARMa Moran's I test for distance-based autocorrelation
## 
## data:  sim
## observed = -0.075229, expected = -0.090909, sd = 0.058432, p-value =
## 0.7884
## alternative hypothesis: Distance-based autocorrelation
# R2 based on Nakagawa et al. 2017: variação explicada pelos efeitos fixos dividida pela total
gab5k.r2 <- r2(gab5k)
gab3k.r2 <- r2(gab3k)

## build AIC table do add model selection results:
aic.tab <- rbind(aic.tab, c("Total abundance FS", 
                            "FC at 5km", 
                            round(gab.fs.aic$AICc[1],2),
                            round(gab.fs.aic$dAICc[1],1), 
                            gab.fs.aic$df[1],
                            round(gab.fs.aic$weight[1],2),
                            round(gab5k.r2$R2_Nagelkerke,2)),
                 c(" ", "~ 1 (NULL)", 
                   round(gab.fs.aic$AICc[2],2), 
                   round(gab.fs.aic$dAICc[2],1),
                   gab.fs.aic$df[2], 
                   round(gab.fs.aic$weight[2],2),
                    "-"),
                 c(" ", "FC at 3km", 
                            round(gab.fs.aic$AICc[3],2),
                            round(gab.fs.aic$dAICc[3],1), 
                            gab.fs.aic$df[3],
                            round(gab.fs.aic$weight[3],2),
                            round(gab3k.r2$R2_Nagelkerke,2)))

Prediction

bRaup-Crick (abundance-based)

Modelling

### LM gamma for regular model selection
rc0 <- lm(brc.abund.fsNA ~ 1, data=data_div)
rc5k <- lm(brc.abund.fsNA ~ fc_5km, data=data_div)
rc5kq <- lm(brc.abund.fsNA ~ fc_5km +  I(fc_5km^2), data=data_div)
rc3k <- lm(brc.abund.fsNA ~ fc_3km, data=data_div)
rc3kq <- lm(brc.abund.fsNA ~ fc_3km +  I(fc_3km^2), data=data_div)

## AIC-based model selection
(rc.fs.aic <- AICctab(rc0, rc5k, rc5kq, 
                      rc3k, rc3kq,
                      base=T,weights=T, logLik=T))
##       logLik AICc dLogLik dAICc df weight
## rc5k   0.1    8.8  3.8     0.0  3  0.595 
## rc3k  -0.8   10.7  2.9     1.8  3  0.237 
## rc0   -3.7   12.7  0.0     3.9  2  0.083 
## rc5kq  0.2   13.3  3.9     4.5  4  0.062 
## rc3kq -0.8   15.4  2.9     6.6  4  0.022

Summary top model:

summary(rc5k)
## 
## Call:
## lm(formula = brc.abund.fsNA ~ fc_5km, data = data_div)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.39189 -0.17995 -0.01784  0.15689  0.50748 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.905410   0.176830   5.120 0.000451 ***
## fc_5km      -0.018211   0.006124  -2.974 0.013955 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2629 on 10 degrees of freedom
## Multiple R-squared:  0.4693, Adjusted R-squared:  0.4162 
## F-statistic: 8.843 on 1 and 10 DF,  p-value: 0.01396
rc2.fs <- rc5k# rename top model to save and plot below

Residuals

sim <- simulateResiduals(rc5k)
plot(sim)

testDispersion(sim) # ok

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.91715, p-value = 0.984
## alternative hypothesis: two.sided
hnp::hnp(rc5k)
## Gaussian model (lm object)

Spatial autocorrelation

testSpatialAutocorrelation(sim, x =  env.data$x, y = env.data$y)

## 
##  DHARMa Moran's I test for distance-based autocorrelation
## 
## data:  sim
## observed = -0.076074, expected = -0.090909, sd = 0.057165, p-value =
## 0.7952
## alternative hypothesis: Distance-based autocorrelation

Prediction

NON-FOREST SPECIALISTS

Gamma diversity

Modelling

### LM alpha for regular model selection
g0 <- glm(gamma_NFS ~ 1, data=data_div, family= Gamma (link = "identity"))
g5k <- glm(gamma_NFS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))
g5kq <- glm(gamma_NFS ~ fc_5km +  I(fc_5km^2), data=data_div, family= Gamma (link = "identity"))
g3k <- glm(gamma_NFS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))
g3kq <- glm(gamma_NFS ~ fc_3km +  I(fc_3km^2), data=data_div, family= Gamma (link = "identity"))


## AIC-based model selection
(g.nfs.aic <- AICctab(g0, 
                      g5k, g5kq, 
                      g3k, g3kq,
                      base=T,weights=T, logLik=T))
##      logLik AICc  dLogLik dAICc df weight
## g3k  -26.3   61.6   7.6     0.0 3  0.7974
## g5k  -28.4   65.9   5.5     4.2 3  0.0958
## g3kq -26.1   65.9   7.8     4.3 4  0.0950
## g5kq -28.4   70.5   5.5     8.9 4  0.0093
## g0   -33.9   73.2   0.0    11.6 2  0.0025

Summary top model:

summary(g3k)
## 
## Call:
## glm(formula = gamma_NFS ~ fc_3km, family = Gamma(link = "identity"), 
##     data = data_div)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.38295  -0.23805  -0.05072   0.13490   0.38893  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 17.39301    1.95020   8.919 4.49e-06 ***
## fc_3km      -0.28151    0.04841  -5.815 0.000169 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.07298721)
## 
##     Null deviance: 2.46666  on 11  degrees of freedom
## Residual deviance: 0.71058  on 10  degrees of freedom
## AIC: 58.646
## 
## Number of Fisher Scoring iterations: 5
g.nfs <- g3k # rename top model to save and plot below

Residuals

sim <- simulateResiduals(g3k)
plot(sim)

testDispersion(sim) # a lil under

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.89592, p-value = 0.936
## alternative hypothesis: two.sided
hnp::hnp(g3k) ## for glm poisson or neg binom (maybe also others)
## Gamma model

boot::glm.diag.plots(g3k) ## for glm poisson or neg binom

Spatial autocorrelation

testSpatialAutocorrelation(sim, x =  env.data$x, y = env.data$y)

## 
##  DHARMa Moran's I test for distance-based autocorrelation
## 
## data:  sim
## observed = -0.123318, expected = -0.090909, sd = 0.057839, p-value =
## 0.5753
## alternative hypothesis: Distance-based autocorrelation

Prediction

Alpha diversity

Modelling

### LM alpha for regular model selection
a0 <- glm(alpha_NFS ~ 1, data=data_div, family= Gamma (link = "identity"))
a5k <- glm(alpha_NFS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))
a5kq <- glm(alpha_NFS ~ fc_5km +  I(fc_5km^2), data=data_div, family= Gamma (link = "identity"))
a3k <- glm(alpha_NFS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))
a3kq <- glm(alpha_NFS ~ fc_3km +  I(fc_3km^2), data=data_div, family= Gamma (link = "identity"))


## AIC-based model selection
(a.nfs.aic <- AICctab(a0, a5k, a5kq, 
                      a3k, a3kq,
                      base=T,weights=T, logLik=T))
##      logLik AICc  dLogLik dAICc df weight
## a3k  -16.1   41.2   5.7     0.0 3  0.497 
## a5k  -16.4   41.8   5.4     0.6 3  0.365 
## a5kq -15.7   45.1   6.1     3.9 4  0.070 
## a3kq -15.9   45.5   5.9     4.3 4  0.057 
## a0   -21.8   49.0   0.0     7.8 2  0.010

Summary top model:

summary(a3k)
## 
## Call:
## glm(formula = alpha_NFS ~ fc_3km, family = Gamma(link = "identity"), 
##     data = data_div)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.66790  -0.14831   0.00666   0.10705   0.54796  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.65928    0.79667   7.104 3.28e-05 ***
## fc_3km      -0.08855    0.02018  -4.388  0.00136 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.1090929)
## 
##     Null deviance: 2.9678  on 11  degrees of freedom
## Residual deviance: 1.1728  on 10  degrees of freedom
## AIC: 38.181
## 
## Number of Fisher Scoring iterations: 5
a.nfs <- a3k # rename top model to save and plot below

Residuals

sim <- simulateResiduals(a3k)
plot(sim)

testDispersion(sim) # a lil under

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.94044, p-value = 0.88
## alternative hypothesis: two.sided
hnp::hnp(a3k) ## for glm poisson or neg binom (maybe also others)
## Gamma model

boot::glm.diag.plots(a3k) ## for glm poisson or neg binom

Spatial autocorrelation

testSpatialAutocorrelation(sim, x =  env.data$x, y = env.data$y)

## 
##  DHARMa Moran's I test for distance-based autocorrelation
## 
## data:  sim
## observed = -0.099144, expected = -0.090909, sd = 0.056241, p-value =
## 0.8836
## alternative hypothesis: Distance-based autocorrelation

Prediction

Beta diversity additive

Modelling

### LM Beta for regular model selection
b0 <- glm(betad_NFS~ 1, data=data_div, family= Gamma (link = "identity"))
b5k <- glm(betad_NFS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))
b5kq <- glm(betad_NFS ~ fc_5km +  I(fc_5km^2), data=data_div, family= Gamma (link = "identity"))
b3k <- glm(betad_NFS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))
b3kq <- glm(betad_NFS ~ fc_3km +  I(fc_3km^2), data=data_div, family= Gamma (link = "identity"))

## AIC-based model selection
(b.nfs.aic <- AICctab(b0, b5k, b5kq, 
                      b3k, b3kq, 
                      base=T,weights=T, logLik=T))
##      logLik AICc  dLogLik dAICc df weight
## b3k  -24.4   57.7   5.5     0.0 3  0.734 
## b3kq -23.8   61.4   6.0     3.7 4  0.118 
## b5k  -26.2   61.4   3.7     3.7 3  0.115 
## b0   -29.9   65.1   0.0     7.4 2  0.019 
## b5kq -26.0   65.6   3.9     7.9 4  0.014

Summary top model:

summary(b3k)
## 
## Call:
## glm(formula = betad_NFS ~ fc_3km, family = Gamma(link = "identity"), 
##     data = data_div)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.71087  -0.18354  -0.04943   0.17316   0.48985  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 11.78046    1.59008   7.409 2.29e-05 ***
## fc_3km      -0.19439    0.03898  -4.987 0.000548 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.1092755)
## 
##     Null deviance: 2.9582  on 11  degrees of freedom
## Residual deviance: 1.2094  on 10  degrees of freedom
## AIC: 54.73
## 
## Number of Fisher Scoring iterations: 6
b.nfs <- b3k # rename top model to save and plot below

Residuals

sim <- simulateResiduals(b3k)
plot(sim)

testDispersion(sim)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.65169, p-value = 0.624
## alternative hypothesis: two.sided
hnp::hnp(b3k) ## for glm poisson or neg binom (maybe also others)
## Gamma model

boot::glm.diag.plots(b3k) ## for glm poisson or neg binom

Spatial autocorrelation

testSpatialAutocorrelation(sim, x =  env.data$x, y = env.data$y)

## 
##  DHARMa Moran's I test for distance-based autocorrelation
## 
## data:  sim
## observed = -0.146700, expected = -0.090909, sd = 0.057505, p-value =
## 0.332
## alternative hypothesis: Distance-based autocorrelation

Prediction

bRaup-Crick (presence/absence-based)

Modelling

### LM gamma for regular model selection
rc0 <- lm(brc.nfsNA ~ 1, data=data_div)
rc5k <- lm(brc.nfsNA ~ fc_5km, data=data_div)
rc5kq <- lm(brc.nfsNA ~ fc_5km +  I(fc_5km^2), data=data_div)
rc3k <- lm(brc.nfsNA ~ fc_3km, data=data_div)
rc3kq <- lm(brc.nfsNA ~ fc_3km +  I(fc_3km^2), data=data_div)

## AIC-based model selection
(rc.nfs.aic <- AICctab(rc0, rc5k, rc5kq, 
                       rc3k, rc3kq,
                       base=T,weights=T, logLik=T))
##       logLik AICc dLogLik dAICc df weight
## rc3k  -0.6   10.1  2.0     0.0  3  0.361 
## rc0   -2.6   10.6  0.0     0.4  2  0.293 
## rc5k  -1.1   11.3  1.5     1.1  3  0.207 
## rc3kq  0.6   12.5  3.2     2.4  4  0.110 
## rc5kq -0.8   15.3  1.8     5.1  4  0.028

Summary of 2nd best model, 1st wqs of absence of effect (~1):

summary(rc3k)
## 
## Call:
## lm(formula = brc.nfsNA ~ fc_3km, data = data_div)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.3020 -0.2156 -0.0678  0.2557  0.4084 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  0.287950   0.190935   1.508   0.1625  
## fc_3km      -0.013092   0.006505  -2.013   0.0719 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2781 on 10 degrees of freedom
## Multiple R-squared:  0.2883, Adjusted R-squared:  0.2171 
## F-statistic: 4.051 on 1 and 10 DF,  p-value: 0.07185
rc1.nfs <- rc3k # rename top model to save and plot below

Residuals

sim <- simulateResiduals(rc3k)
plot(sim)

testDispersion(sim) # a lil under

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.91715, p-value = 0.984
## alternative hypothesis: two.sided
hnp::hnp(rc3k)
## Gaussian model (lm object)

Spatial autocorrelation

testSpatialAutocorrelation(sim, x =  env.data$x, y = env.data$y)

## 
##  DHARMa Moran's I test for distance-based autocorrelation
## 
## data:  sim
## observed = -0.159951, expected = -0.090909, sd = 0.058625, p-value =
## 0.2389
## alternative hypothesis: Distance-based autocorrelation

Prediction

Total abundance

Modelling

### LM "gamma"  abundanec for regular model selection
gab0 <- glm.nb(ab.land_NFS ~ 1, data=data_div)
gab5k <- glm.nb(ab.land_NFS ~ fc_5km, data=data_div)
gab5kq <- glm.nb(ab.land_NFS ~ fc_5km +  I(fc_5km^2), data=data_div)
gab3k <- glm.nb(ab.land_NFS ~ fc_3km, data=data_div)
gab3kq <- glm.nb(ab.land_NFS ~ fc_3km +  I(fc_3km^2), data=data_div)

## AIC-based model selection
(gab.nfs.aic <- AICctab(gab0, gab5k, gab5kq, 
                        gab3k, gab3kq,
                        base=T,weights=T, logLik=T))
##        logLik AICc  dLogLik dAICc df weight
## gab5k  -63.9  136.7   5.3     0.0 3  0.402 
## gab3k  -64.2  137.4   4.9     0.7 3  0.285 
## gab5kq -62.0  137.7   7.1     0.9 4  0.253 
## gab3kq -63.6  141.0   5.5     4.3 4  0.047 
## gab0   -69.1  143.6   0.0     6.8 2  0.013

Summary top model:

summary(gab5k)
## 
## Call:
## glm.nb(formula = ab.land_NFS ~ fc_5km, data = data_div, init.theta = 4.180831522, 
##     link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.8524  -0.4383   0.2265   0.4858   0.8634  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  5.92551    0.33598   17.64  < 2e-16 ***
## fc_5km      -0.04812    0.01180   -4.08  4.5e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(4.1808) family taken to be 1)
## 
##     Null deviance: 29.171  on 11  degrees of freedom
## Residual deviance: 12.765  on 10  degrees of freedom
## AIC: 133.73
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  4.18 
##           Std. Err.:  1.78 
## 
##  2 x log-likelihood:  -127.727
gab.nfs <- gab5k # rename top model to save and plot below

Residuals

sim <- simulateResiduals(gab5k)
plot(sim)

testDispersion(sim) # a lil under

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.43345, p-value = 0.192
## alternative hypothesis: two.sided
hnp::hnp(gab5k)
## Negative binomial model (using MASS package)

boot::glm.diag.plots(gab5k)

Spatial autocorrelation

testSpatialAutocorrelation(sim, x =  env.data$x, y = env.data$y)

## 
##  DHARMa Moran's I test for distance-based autocorrelation
## 
## data:  sim
## observed = -0.079003, expected = -0.090909, sd = 0.055904, p-value =
## 0.8313
## alternative hypothesis: Distance-based autocorrelation
# R2 based on Nakagawa et al. 2017: variação explicada pelos efeitos fixos dividida pela total
gab5k.r2 <- r2(gab5k)
gab5kq.r2 <- r2(gab5kq)
gab3k.r2 <- r2(gab3k)

## build AIC table do add model selection results:
aic.tab <- rbind(aic.tab, c("Total abundance NFS", 
                            "FC at 5km", 
                            round(gab.nfs.aic$AICc[1],2),
                            round(gab.nfs.aic$dAICc[1],1), 
                            gab.nfs.aic$df[1],
                            round(gab.nfs.aic$weight[1],2),
                            round(gab5k.r2$R2_Nagelkerke,2)),
                 c(" ", "FC at 3km",
                   round(gab.nfs.aic$AICc[2],2), 
                   round(gab.nfs.aic$dAICc[2],1),
                   gab.nfs.aic$df[2], 
                   round(gab.nfs.aic$weight[2],2),
                   round(gab3k.r2$R2_Nagelkerke,2)),
                 c(" ", "Non-linear FC at 5km",
                   round(gab.nfs.aic$AICc[3],2), 
                   round(gab.nfs.aic$dAICc[3],1),
                   gab.nfs.aic$df[3], 
                   round(gab.nfs.aic$weight[3],2),
                   round(gab5kq.r2$R2_Nagelkerke,2)
                 ))

Prediction

bRaup-Crick (abundance-based)

Modelling

### LM gamma for regular model selection
rc0 <- lm(brc.abund.nfsNA ~ 1, data=data_div)
rc5k <- lm(brc.abund.nfsNA ~ fc_5km, data=data_div)
rc5kq <- lm(brc.abund.nfsNA ~ fc_5km +  I(fc_5km^2), data=data_div)
rc3k <- lm(brc.abund.nfsNA ~ fc_3km, data=data_div)
rc3kq <- lm(brc.abund.nfsNA ~ fc_3km +  I(fc_3km^2), data=data_div)

## AIC-based model selection
(rc.nfs.aic <- AICctab(rc0, rc5k, rc5kq, 
                       rc3k, rc3kq,
                       base=T,weights=T, logLik=T))
##       logLik AICc dLogLik dAICc df weight
## rc3k  -7.4   23.8  1.9     0.0  3  0.367 
## rc0   -9.3   23.9  0.0     0.1  2  0.356 
## rc5k  -8.0   24.9  1.3     1.1  3  0.215 
## rc3kq -7.3   28.3  2.0     4.5  4  0.040 
## rc5kq -7.9   29.5  1.4     5.6  4  0.022

Summary 2nd best model, 1st is of absence of effect (~1):

summary(rc3k)
## 
## Call:
## lm(formula = brc.abund.nfsNA ~ fc_3km, data = data_div)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.23916 -0.18519  0.05082  0.27801  0.53338 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  0.82320    0.33770   2.438   0.0350 *
## fc_3km      -0.02195    0.01151  -1.908   0.0855 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4918 on 10 degrees of freedom
## Multiple R-squared:  0.2669, Adjusted R-squared:  0.1936 
## F-statistic: 3.641 on 1 and 10 DF,  p-value: 0.08547
rc2.nfs <- rc3k # rename top model to save and plot below

Residuals

sim <- simulateResiduals(rc3k)
plot(sim)
## qu = 0.75, log(sigma) = -3.631427 : outer Newton did not converge fully.

testDispersion(sim) # a lil under

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.91715, p-value = 0.984
## alternative hypothesis: two.sided
hnp::hnp(rc3k)
## Gaussian model (lm object)

Spatial autocorrelation

testSpatialAutocorrelation(sim, x =  env.data$x, y = env.data$y)

## 
##  DHARMa Moran's I test for distance-based autocorrelation
## 
## data:  sim
## observed = -0.118475, expected = -0.090909, sd = 0.055483, p-value =
## 0.6193
## alternative hypothesis: Distance-based autocorrelation

Prediction

AICc-model selection table

Table S2. AICc-based selected models table containing the models of gamma (γ), alpha (ɑ), additive beta (β) diversities followed by abundance at landscape total, abundance at mean site and β-diversity calculated from the Raup-Crick null model approach, for both presence/absence (βRC) and abundance-based (βRC-abundance) that scored below ΔAICc=2. All tested as a function of Atlantic Forest loss for each group of species, first models separated into forest specialists species (FS; Table 1a), and followed by non-forest species (NFS; Table 1b). We show selected models ordered from the lowest to the highest AICc values ΔAICc≤ 2. We present information regarding the diversity metric tested (Diversity), the model terms for forest loss predictors (Model), the Akaike information criterion corrected for small samples (AICc), the AICc difference from the first-ranked model (ΔAICc), degrees of freedom (df) Akaike weight (ωi) and model variance explained (r²).
Diversity Model AICc dAICc df weight r2
Gamma FS FC at 5km 65.36 0 3 0.75 0.58
Alpha FS FC at 5km 47.99 0 3 0.54 0.4
FC at 3km 49.97 2 3 0.2 0.29
Beta FS FC at 5km 53.97 0 3 0.71 0.61
RCbeta (P/A-based) FS ~ 1 (NULL) -4.25 0 2 0.42
Non-linear FC at 3km -3.02 1.2 4 0.23 0.33
Total abundance FS FC at 5km 150.13 0 3 0.4 0.44
~ 1 (NULL) 150.66 0.5 2 0.31
FC at 3km 151.37 1.2 3 0.22 0.33
RCbeta (abund-based) FS FC at 5km 8.81 0 3 0.6 0.42
FC at 3km 10.65 1.8 3 0.24 0.32
Gamma NFS FC at 3km 61.65 0 3 0.8 0.73
Alpha NFS FC at 3km 41.18 0 3 0.5 0.63
FC at 5km 41.8 0.6 3 0.37 0.61
Beta NFS FC at 3km 57.73 0 3 0.73 0.62
RCbeta (P/A-based) NFS FC at 3km 10.15 0 3 0.36 0.22
~ 1 (NULL) 10.56 0.4 2 0.29
FC at 5km 11.26 1.1 3 0.21 0.14
Total abundance NFS FC at 5km 136.73 0 3 0.4 0.82
FC at 3km 137.42 0.7 3 0.28 0.79
Non-linear FC at 5km 137.66 0.9 4 0.25 0.93
RCbeta (abund-based) NFS FC at 3km 23.83 0 3 0.37 0.19
~ 1 (NULL) 23.89 0.1 2 0.36
FC at 5km 24.9 1.07 3 0.22 0.12

MODEL COEFFICIENTS TABLE

# fazer tabela de coeficientes dos modelos selectionados
(tab_gfs <- tab_model(g.fs, show.p = F, show.se = TRUE))
  gamma FS
Predictors Estimates std. Error CI
(Intercept) 21.47 2.02 17.69 – 25.70
fc 5km -0.23 0.06 -0.35 – -0.11
Observations 12
R2 Nagelkerke 0.579 
(tab_afs <- tab_model(a.fs, show.p = F, show.se = TRUE))
  alpha FS
Predictors Estimates std. Error CI
(Intercept) 7.57 0.96 5.84 – 9.57
fc 5km -0.08 0.03 -0.13 – -0.02
Observations 12
R2 Nagelkerke 0.396 
(tab_bfs <- tab_model(b.fs, show.p = F, show.se = TRUE))
  betad FS
Predictors Estimates std. Error CI
(Intercept) 13.94 1.28 11.50 – 16.66
fc 5km -0.16 0.04 -0.23 – -0.08
Observations 12
R2 Nagelkerke 0.613 
(tab_brc1fs <- tab_model(rc1.fs, show.p = F, show.se = TRUE))
  brc fs NA
Predictors Estimates std. Error CI
(Intercept) -0.23 0.05 -0.33 – -0.12
Observations 12
R2 / R2 adjusted 0.000 / 0.000 
(tab_gabfs <- tab_model(gab.fs, show.p = F, show.se = TRUE))
  ab land FS
Predictors Incidence Rate Ratios std. Error CI
(Intercept) 379.43 120.98 208.54 – 714.06
fc 5km 0.98 0.01 0.96 – 1.00
Observations 12
R2 Nagelkerke 0.442 
(tab_brc2fs <- tab_model(rc2.fs, show.p = F, show.se = TRUE))
  brc abund fs NA
Predictors Estimates std. Error CI
(Intercept) 0.91 0.18 0.51 – 1.30
fc 5km -0.02 0.01 -0.03 – -0.00
Observations 12
R2 / R2 adjusted 0.469 / 0.416 
    # Nnfs
(tab_gnfs <- tab_model(g.nfs, show.p = F, show.se = TRUE))
  gamma NFS
Predictors Estimates std. Error CI
(Intercept) 17.39 1.95 13.79 – 21.65
fc 3km -0.28 0.05 -0.38 – -0.18
Observations 12
R2 Nagelkerke 0.733 
(tab_anfs <- tab_model(a.nfs, show.p = F, show.se = TRUE))
  alpha NFS
Predictors Estimates std. Error CI
(Intercept) 5.66 0.80 4.28 – 7.39
fc 3km -0.09 0.02 -0.13 – -0.05
Observations 12
R2 Nagelkerke 0.634 
(tab_bnfs <- tab_model(b.nfs, show.p = F, show.se = TRUE))
  betad NFS
Predictors Estimates std. Error CI
(Intercept) 11.78 1.59 8.81 – 15.36
fc 3km -0.19 0.04 -0.28 – -0.11
Observations 12
R2 Nagelkerke 0.621 
(tab_brc1nfs <- tab_model(rc1.nfs, show.p = F, show.se = TRUE))
  brc nfs NA
Predictors Estimates std. Error CI
(Intercept) 0.29 0.19 -0.14 – 0.71
fc 3km -0.01 0.01 -0.03 – 0.00
Observations 12
R2 / R2 adjusted 0.288 / 0.217 
(tab_gabnfs <- tab_model(gab.nfs, show.p = F, show.se = TRUE))
  ab land NFS
Predictors Incidence Rate Ratios std. Error CI
(Intercept) 374.47 125.81 204.76 – 711.73
fc 5km 0.95 0.01 0.93 – 0.97
Observations 12
R2 Nagelkerke 0.817 
(tab_brc2nfs <- tab_model(rc2.nfs, show.p = F, show.se = TRUE))
  brc abund nfs NA
Predictors Estimates std. Error CI
(Intercept) 0.82 0.34 0.07 – 1.58
fc 3km -0.02 0.01 -0.05 – 0.00
Observations 12
R2 / R2 adjusted 0.267 / 0.194

FIGURES

Gamma hab

### Gamma ####
# r2(g.fs) # 0.55
## prediction g2
new <- data.frame(fc_5km = seq(min(data_div$fc_5km),max(data_div$fc_5km), length.out=50))
preds <- predict(g.fs, type="response", newdata=new, se.fit=T)

new$pred <- preds$fit
new$ic.up <- preds$fit + preds$se.fit*1.96 
new$ic.low <- preds$fit - preds$se.fit*1.96 

# r2(g.nfs) # 0.73

## prediction g1
new2 <- data.frame(fc_3km = seq(min(data_div$fc_3km),max(data_div$fc_3km), length.out=50))
preds2 <- predict(g.nfs, type="response", newdata=new2, se.fit=T)

new2$pred <- preds2$fit
new2$ic.up <- preds2$fit + preds2$se.fit*1.96 
new2$ic.low <- preds2$fit - preds2$se.fit*1.96 

## Plot gamma richness ~ fc %
(g.hab <- ggplot(data= new, aes(x= fc_5km, y = pred)) + 
    geom_line(data= new, aes(x=fc_5km, y=pred), col= "#0072B2") +
    geom_ribbon(data= new, aes(x=fc_5km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#0072B2") +
    geom_point(data= data_div, aes(x=fc_5km, y=gamma_FS), col= "#0072B2") + 
    geom_ribbon(data= new2, aes(x=fc_3km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#D55E00") +
    geom_line(data= new2, aes(x=fc_3km, y=pred), col= "#D55E00") +
    geom_point(data= data_div, aes(x=fc_3km, y=gamma_NFS), col= "#D55E00") +
    xlab(" ") + ylab("Gamma") + 
    theme(axis.text=element_text(size=15),axis.title=element_text(size=14,face="bold"),legend.position = 'none') + 
    ylim(0,max(data_div[2:7])+0.5) + 
    # scale_x_reverse() +
    annotate("text", y=20.5, x=38,colour= "#0072B2",
             label=expression(paste(~ R^2 ~ "= 0.55"))) +
    annotate("text", y=19, x=40, colour= "#D55E00",
             label=expression(paste(~ R^2 ~ "= 0.73"))) + labs(tag="a"))

Alpha hab

# r2(a.fs) # 0.40

## prediction a2
new <- data.frame(fc_5km = seq(min(data_div$fc_5km),max(data_div$fc_5km), length.out=50))
preds <- predict(a.fs, type="response", newdata=new, se.fit=T)

new$pred <- preds$fit
new$ic.up <- preds$fit + preds$se.fit*1.96 
new$ic.low <- preds$fit - preds$se.fit*1.96 

# r2(a.nfs) # 0.63
## prediction a1
new2 <-  data.frame(fc_3km = seq(min(data_div$fc_3km),max(data_div$fc_3km), length.out=50))
preds <- predict(a.nfs, type="response", newdata=new2, se.fit=T)

new2$pred <- preds$fit
new2$ic.up <- preds$fit + preds$se.fit*1.96 
new2$ic.low <- preds$fit - preds$se.fit*1.96 

## Plot Alpha richness ~ fc_3km
(a.hab <- ggplot(data= new, aes(x=fc_5km, y=pred)) + 
  geom_line(data= new, aes(x=fc_5km, y=pred), col= "#0072B2") +
  geom_ribbon(data= new, aes(x=fc_5km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#0072B2") +
  geom_point(data=data_div, aes(x=fc_5km, y=alpha_FS), col= "#0072B2") +
  geom_line(data= new2, aes(x=fc_3km, y=pred), col= "#D55E00") +
  geom_ribbon(data= new2, aes(x=fc_3km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#D55E00") +
  geom_point(data=data_div, aes(x=fc_3km, y=alpha_NFS), col= "#D55E00") +
  xlab("") + ylab("Alpha") + 
  scale_color_manual(breaks=c("FS", "NFS"),
                     values=c('FS'='#0072B2', 'NFS'='#D55E00')) +
  theme(axis.text=element_text(size=15),
        axis.title=element_text(size=14,face="bold"),
        legend.position = 'bottom') + 
  ylim(0,max(data_div[2:7])+0.5) +
    # scale_x_reverse() +
  annotate("text", y=13, x=38,colour= "#0072B2",
           label=expression(paste(~ R^2 ~ "= 0.40")))+
  annotate("text", y=11, x=40, colour= "#D55E00",
           label=expression(paste(~ R^2 ~ "= 0.63"))) +
  labs(tag="b"))

Beta hab

### Beta add  ####
# r2(b.fs) # 0.58

## prediction b1
new <- data.frame(fc_5km = seq(min(data_div$fc_5km),max(data_div$fc_5km), length.out=50))
preds <- predict(b.fs, type="response", newdata=new, se.fit=T)

new$pred <- preds$fit
new$ic.up <- preds$fit + preds$se.fit*1.96 
new$ic.low <- preds$fit - preds$se.fit*1.96 

# r2(b.nfs) # 0.65

## prediction b1
new2 <-  data.frame(fc_3km = seq(min(data_div$fc_3km),max(data_div$fc_3km), length.out=50))
preds <- predict(b.nfs, type="response", newdata=new2, se.fit=T)

new2$pred <- preds$fit
new2$ic.up <- preds$fit + preds$se.fit*1.96 
new2$ic.low <- preds$fit - preds$se.fit*1.96 

## Plot Beta Diversity ~ fc %
(b.hab <- ggplot(data= new, aes(x=fc_5km, y=pred)) + 
  geom_line(data= new, aes(x=fc_5km, y=pred), col= "#0072B2") +
  geom_ribbon(data= new, aes(x=fc_5km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#0072B2") +
  geom_point(data=data_div, aes(x=fc_5km, y=betad_FS), col= "#0072B2") +
  geom_line(data= new2, aes(x= fc_3km, y= pred), 
            col= "#D55E00") +
  geom_ribbon(data= new2, aes(x=fc_3km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#D55E00") +
  geom_point(data=data_div, aes(x=fc_3km, y=betad_NFS), col= "#D55E00") +
  xlab(" ") + ylab("Beta") + 
  theme(axis.text=element_text(size=15),axis.title=element_text(size=14,face="bold"),legend.position = 'none') + 
  ylim(0,max(data_div[2:7])+0.5) + 
    # scale_x_reverse() +
  annotate("text", y=18, x=40,colour= "#0072B2",
           label=expression(paste(~ R^2 ~ "= 0.58")))+
  annotate("text", y=16, x=38, colour= "#D55E00",
           label=expression(paste(~ R^2 ~ "= 0.65"))) +
  labs(tag="c"))

g.hab + a.hab + b.hab

ggsave("figures/plot.hab.jpeg",units="cm", width=20, height=8, dpi=300)

Total abundance FS & NFS

## FS
# r2(gab.fs) # 0.44

## prediction gab1
new <- data.frame(fc_5km = seq(min(data_div$fc_5km),max(data_div$fc_5km), length.out=50))
preds <- predict(gab.fs, type="response", newdata=new, se.fit=T)

new$pred <- preds$fit
new$ic.up <- preds$fit + preds$se.fit*1.96 
new$ic.low <- preds$fit - preds$se.fit*1.96 

## NFS
# r2(gab.nfs) # 0.82

## prediction gab2
new2 <- new
preds <- predict(gab.nfs, type="response", newdata=new2, se.fit=T)

new2$pred <- preds$fit
new2$ic.up <- preds$fit + preds$se.fit*1.96 
new2$ic.low <- preds$fit - preds$se.fit*1.96 

## Plot gamma richness ~ forest cover 5km
(gab.hab <- ggplot(new, aes(x=fc_5km, y=pred)) + 
  geom_line(linetype= "dashed", col= "#0072B2", aes(alpha= 0.1)) +
  # geom_ribbon(aes(x=fc_5km , ymax=ic.up, ymin=ic.low), alpha=0.1,  fill= "#0072B2") +
  geom_point(data=data_div, aes(x=fc_5km, y=ab.land_FS), col= "#0072B2") +
  geom_line(data= new2, aes(x=fc_5km, y=pred), col= "#D55E00") +
  geom_ribbon(data= new2, aes(x=fc_5km , ymax=ic.up, ymin=ic.low), alpha=0.2, fill= "#D55E00") +
  geom_point(data= data_div, aes(x=fc_5km, y=ab.land_NFS), col= "#D55E00") +
    scale_x_continuous(limits= c(10,50)) +
  xlab(" ") + ylab("Total abundance") + 
  theme(axis.text=element_text(size=15),axis.title=element_text(size=14,face="bold"),legend.position = 'none') + 
  # scale_x_reverse() +
    annotate("text", y=300, x=38, alpha= .5, colour= "#0072B2",
             label=expression(paste(~ R^2 ~ "= 0.44"))) +
    annotate("text", y=270, x=40, colour= "#D55E00",
           label=expression(paste(~ R^2 ~ "= 0.82"))) +
  labs(tag="d"))

### Beta RC FS & NFS

## FS
# ## NULL ## r2(rc1.fs) 

## prediction rc2.2 and rc2
new <- data.frame(fc_3km = seq(min(data_div$fc_3km),max(data_div$fc_3km), length.out=30))
preds <- predict(rc1.fs, type="response", newdata=new, se.fit=T)

new$pred <- preds$fit
new$ic.up <- preds$fit + preds$se.fit*1.96 
new$ic.low <- preds$fit - preds$se.fit*1.96 

## NFS
# r2(rc1.nfs) # 0.25

## prediction rc2
new2 <- data.frame(fc_3km = seq(min(data_div$fc_3km),max(data_div$fc_3km), length.out=50))
preds <- predict(rc1.nfs, type="response", newdata=new2, se.fit=T)

new2$pred <- preds$fit
new2$ic.up <- preds$fit + preds$se.fit*1.96 
new2$ic.low <- preds$fit - preds$se.fit*1.96 

## Plot RCbeta abundance ~ forest cover at 3
(rc.hab <- ggplot(data= new, aes(x=fc_3km, y=pred)) + 
  geom_line(data= new,aes(x=fc_3km, y=pred), 
            col= "#0072B2", linetype= "dashed", alpha= .2) +
  # geom_ribbon(data= new, aes(x=fc_3km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#0072B2") + 
  geom_point(data=data_div, aes(x=fc_3km, y=brc.fsNA), col= "#0072B2") +
  geom_line(data= new2, aes(x=fc_3km, y=pred), col= "#D33E00", linetype= "dashed", alpha= .2) +
  # geom_ribbon(data= new2, aes(x=fc_3km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#D33E00") +
  geom_point(data=data_div, aes(x=fc_3km, y=brc.nfsNA), col= "#D33E00") +
  xlab("Percent forest cover") + ylab(c(expression(bold("β"["RC"])))) + 
  theme(axis.text=element_text(size=14), axis.title=element_text(size=14,face="bold"),
        legend.position = 'none') + ylim(-1,1) + 
    # scale_x_reverse() +
  annotate("text", y=0.9, x=43, alpha= .5, colour= "#D33E00",
           label=expression(paste(~ R^2 ~ "= 0.25"))) +
  # annotate("text", y=0.7, x=43, alpha= .5, colour= "#0072B2",
           # label=expression(paste(~ R^2 ~ "= 0.30"))) +
  # annotate("text", y=0.9, x=43, colour= "#D33E00",
  #          label=expression(paste(beta, "(SE) = -0.06(0.03); " ~ R^2 ~ "= 0.37"))) +
  # annotate("text", y=0.8, x=43, colour= "#0072B2",
  #          label=expression(paste(beta, "(SE) = -0.06(0.03); " ~ R^2 ~ "= 0.33"))) +
  labs(tag="e"))

Beta RC abund FS & NFS

# r2(rc2.fs) # 0.41

## prediction rc2
new <- data.frame(fc_5km = seq(min(data_div$fc_5km),max(data_div$fc_5km), length.out=50))
preds <- predict(rc2.fs, type="response", newdata=new, se.fit=T)

new$pred <- preds$fit
new$ic.up <- preds$fit + preds$se.fit*1.96 
new$ic.low <- preds$fit - preds$se.fit*1.96 

## NFS
# r2(rc2.nfs) # 0.21

## prediction rc1
new2 <- data.frame(fc_3km = seq(min(data_div$fc_3km),max(data_div$fc_3km), length.out=50))
preds <- predict(rc2.nfs, type="response", newdata=new2, se.fit=T)

new2$pred <- preds$fit
new2$ic.up <- preds$fit + preds$se.fit*1.96 
new2$ic.low <- preds$fit - preds$se.fit*1.96 

## RCbeta abund ~ FC%
(rc.ab.hab <- ggplot(data= new, aes(x=fc_5km, y=pred)) + 
  geom_line(data= new,aes(x=fc_5km, y=pred), col= "#0072B2") +
  geom_ribbon(data= new, aes(x=fc_5km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#0072B2") + 
  geom_point(data=data_div, aes(x=fc_5km, y=brc.abund.fsNA), col= "#0072B2") +
  geom_line(data= new2, aes(x=fc_3km, y=pred), linetype= "dashed", col= "#D55E00", alpha=.2) +
  # geom_ribbon(data= new2, aes(x=fc_3km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#D55E00") +
  geom_point(data=data_div, aes(x=fc_3km, y=brc.abund.nfsNA), col= "#D55E00") +
  xlab(" ") + ylab(c(expression(bold("β"["RC-abundance"])))) + 
  theme(axis.text=element_text(size=14), axis.title=element_text(size=14,face="bold"),legend.position = 'none') +
    ylim(-1,1.1) + 
    # scale_x_reverse() +
  annotate("text", y=0.9, x=45, colour= "#0072B2",
           label=expression(paste(~ R^2 ~ "= 0.41"))) +
  annotate("text", y=0.75, x=43, alpha= .5, colour= "#D55E00",
           label=expression(paste( ~ R^2 ~ "= 0.21"))) +
  labs(tag="f"))

Final plot

(plot.hab <- (g.hab + a.hab + b.hab) /
    (gab.hab + rc.hab + rc.ab.hab))

ggsave("figures/figure2.jpeg",units="cm", width=28, height=15, dpi=300)